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Abstract. We investigate the classical counterpart of an effective Hamiltonian for a 
strongly trimerized kagome lattice. Ahhough the Hamiltonian only has a discrete 
symmetry, the classical groundstate manifold has a continuous global rotational 
symmetry. Two cases should be distinguished for the sign of the exchange constant. 
In one case, the groundstate has a 120° spin structure. To determine the transition 
temperature, we perform Monte-Carlo simulations and measure specific heat, the order 
parameter as well as the associated Binder cumulant. In the other case, the classical 
groundstates are macroscopically degenerate. A thermal order-by-disorder mechanism 
is predicted to select another 120° spin-structure. A finite but very small transition 
temperature is detected by Monte-Carlo simulations using the exchange method. 
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1. Introduction 

Highly frustrated magnets are a fascinating area of research with many challenges and 
surprises [1,2]. One exotic case is the spin- 1/2 Heisenberg antiferromagnet on the 
kagome lattice, where in particular an unusually large number of low-lying singlets 
was observed numerically [3,4]. An effective Hamiltonian approach to the strongly 
trimerized kagome lattice [5] was then successful in explaining the unusual properties of 
the homogeneous lattice [6]. 

Interest in this situation has been renewed recently due to the suggestion that 
strongly trimerized kagome lattices can be realized by fermionic quantum gases in 
optical lattices [7,8]. For two fermions per triangle, one obtains the aforementioned 
effective Hamiltonian on a triangular lattice [5], but without the original magnetic 
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Figure 1. Part of the triangular lattice. Arrows indicate the direction of the unit 
vectors ei-{ij) entering the Hamiltonian 

degrees of freedom. This effective Hamiltonian also describes the spin-1/2 Heisenberg 
antiferromagnet on the trimerized kagome lattice at one third of the saturation 
magnetization [9]. Furthermore, this model shares some features with pure orbital 
models on the square lattice [10,11], but it is substantially more frustrated than these 
models. 

Numerical studies of the effective quantum Hamiltonian [12-14] provide evidence for 
an ordered groundstate and a possible finite-temperature ordering transition. Since the 
Hamiltonian only has discrete symmetries, one expects indeed a finite-temperature phase 
transition if the groundstate is ordered. Furthermore, quantum fluctuations should be 
unimportant for the generic properties of such a phase transition, which motivates us 
to study the classical counterpart of the model at finite temperatures. 

2. Model and symmetries 

In this paper we study a Hamiltonian which is given in terms of spins Si by 

^='^J2{^ ^i'i^'i) ■ (2 ^r,{^,l) ■ Sj) ■ (1) 

The sum runs over the bonds of nearest neighbours {i,j) of a triangular lattice with 
N sites. For each bond, only certain projections of the spins Si on unit vectors ei-{ij) 
enter the interaction. These directions are sketched in Fig. ^ Note that these directions 
depend both on the bond {i,j) and the corresponding end i or j such that three different 
projections of each spin Si enter the interaction with its six nearest neighbours. 

In the derivation from the trimerized kagome lattice [5,7, 15], the Si are pseudo- 
spin operators acting on the two chiralities on each triangle and should therefore be 
considered as quantum spin-1/2 operators. Here we will treat the Si as classical unit 
vectors. Since only the x- and y-components enter the Hamiltonian (^, one may take 
the Si as 'planar' two-component vectors. On the other hand, in the quantum case 
commutation relations dictate the presence of the 2;-component as well, such that taking 
the Si as 'spherical' three-component vectors is another natural choice [14]. Here we 
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will compare both choices and thus assess qualitatively the effect of omitting the z- 
components. 

The internal symmetries of the Hamiltonian (0) constitute the dihedral group Dq, 
i.e., the symmetry group of a regular hexagon. Some of its elements consist of a 
simultaneous transformation of the spins Si and the lattice [14]. Since these symmetries 
are only discrete, a finite-temperature phase transition is allowed above an ordered 
groundstate. 

The derivation from a spin model [5, 15] yields a positive J > 0, while it may also 
be possible to realize J < [8] in a Fermi gas in an optical lattice [7]. In the following 
we will first discuss the case of a negative exchange constant J < and then turn to 
the case of a positive exchange constant J > 0. The second case is more interesting, 
but also turns out to be more difficult to handle. 

3. Negative exchange constant 

For J < there is a one-parameter family of ordered groundstates. Si = 
(cos(0j), sin(0j), 0) with (pa = 0, 4>b = + 2it/3 and 0c = — 27i/3 on the three 
sublattices a, b, and c, respectively (120° Neel order). The energy of these states, 
£]J<o — Qj]\f^ ig independent of 6. Computing the free energy JF"^<°(6') by including 
the effect of Gaussian fluctuations we find that J^'^^^{6) has minima at ^ = (2 n-|- 1) vr/G, 
n = 0, 1, . . . , 5 [14]. This implies that the above 120° Neel structures lock in at these 
angles. 

We have performed Monte-Carlo (MC) simulations for J < using a standard 
single-spin flip Metropolis algorithm [16]. A snapshot of a low-temperature configuration 
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Figure 3. MC results for the specific heat C for J < 0. FuU lines are for planar spins; 
dashed lines for spherical spins. Increasing line widths denote increasing system sizes 
N. Error bars are negligible on the scale of the figure. 



on an = 12 X 12 lattice is shown in Fig. |2ta). The 120° ordering is clearly seen in 
such snapshots. Fig. Efb) shows histograms of the angles 0j. One observes that with 
increasing lattice size N, pronounced maxima emerge in the probability P{(f)i) to observe 
an angle (pi at the predicted lock-in values 6 = (2n + 1) 7r/6. 

Thermodynamic quantities have been computed by averaging over at least 100 
independent MC simulations. Each simulation was started at high temperatures and 
slowly cooled to lower temperatures in order to minimize equilibration times. A first 
quantity, namely the specific heat C, is shown in Fig. 01 There is a maximum in C at 
T ^ 1.8 \J\ for spherical spins, and for planar spins at a higher temperature T ^ 2.2 \J\. 
The fact that the value of C around the maximum increases with N indicates a phase 
transition. For T ^ 0, the equipartition theorem predicts a contribution 1/2 to the 
specific heat per transverse degree of freedom. Indeed, the low-temperature results are 
very close to C/N = 1/2 and 1 for planar and spherical spins, respectively. 

In order to quantify the expected order, we introduce the sublattice order parameter 



where the sum runs over one of the three sublattices C of the triangular lattice. Fig. 0] 
plots the square of this sublattice order parameter 



which is a scalar quantity and expected to be non-zero in a three-sublattice ordered state. 
Note, however, that the order parameter (0) is insensitive to the mutual orientation of 
spins on the three sublattices. In Fig. 0] we observe that converges to zero with 
oo at high temperatures, while a non-zero value persists at lower temperatures. 




(2) 




(3) 




Figure 4. MC results for the square of the sublattice magnetization for J < 0. 
Full lines are for planar spins; dashed lines for spherical spins. Increasing line widths 
denote increasing N. Error bars are negligible on the scale of the figure. 



consistent with a phase transition around T ^ 2 | J|, as aheady indicated by the specific 
heat. Furthermore, the sublattice order again points to a higher transition temperature 
for planar spins than for spherical spins. It should also be noted that there are noticeable 
quantitative differences between the values of for planar and spherical spins over the 
entire temperature range. Just for T ^ both planar and spherical spins yield — 1, 
as expected for a perfectly ordered groundstate. 

In order to determine the transition temperature more accurately, we use the 
'Binder cumulant' [16, 17] associated to the order parameter (j2)) via 



( ) f 1 for planar spins. 



U4 = l + A-A ^ % with A = < 2, (4) 

l^M'i^ y - for spherical spins. 

The constants in (0)) are chosen such that = for a Gaussian distribution of the 
order parameter P{Ms) oc exp ^— cM^^ Such a distribution is expected at high 
temperatures, leading to f/4 — > for T ^ \J\. Conversely, a perfectly ordered state 
yields {Mf) = (Mf / such that with the prefactors as in (jH) we find U4 = 1. Hence, 



for an ordered state we expect ?74 ~ 1 for T < Tc. 

Fig. El shows our results for the Binder cumulant, as defined in (jl)). The transition 
temperature can be estimated from the crossings of the Binder cumulants at different 
sizes [16, 17], which are shown for planar spins in the inset of Fig. |S] Since we 
have not been aiming at high precision, we cannot perform a finite-size extrapolation. 
Nevertheless, we obtain a rough estimate 

yplanar^2|J|, (5) 

with an error on the order of a few percent. The corresponding value for spherical 
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Figure 5. Main panel: MC results for the Binder cumulant for J < 0. Full lines are for 
planar spins; dashed lines for spherical spins. Increasing line widths denote increasing 
N. Error bars are negligible on the scale of the figure. Inset: Binder cumulant for 
planar spins close to the critical temperature. 

spins is estimated as T^piif^^^cai ^ 1.57 |J| [14]. We therefore conclude that out-of-plane 
fluctuations reduce Tc by about 20%. 

4. Positive exchange constant 

As in the case J < 0, there is a one-parameter family of 120° Neel ordered groundstates 
with energy —3 J N. They differ from the states found for J < by an interchange 
of the spin directions on the b- and c-sublattices. The contribution of Gaussian 
fluctuations around these states yields a free energy J^^^^{9) which has minima at 
9 = nn/S, n = 0, 1, . . . , 5. Hence the Neel structure locks in at these angles for 
J > 0. Surprisingly, an inspection of all states of finite cells of the lattice, in which the 
mutual angles between pairs of spins are multiples of 2tc/3 [14], reveals that there is a 
macroscopic number of groundstates in addition to the 120° Neel states, i.e., the number 
of ground states grows exponentially with A^. The 120° Neel state has soft modes, 
all other groundstates have fewer soft modes [14]. At low but finite temperatures the 
120° Neel state will therefore be selected by a thermal order- by-disorder mechanism [18]. 

The large number of groundstates which are separated by an energy barrier 
renders it extremely difficult to thermalize a simple MC simulation at sufficiently low 
temperatures. We have therefore performed exchange MC simulations (also known 
as 'parallel tempering') [19, 20] for J > 0, using a parallel implementation based on 
MPI. 96 replicas were distributed over the temperature range T/J = 0.01, . . . , 0.7 in a 
manner to ensure an acceptance rate for exchange moves of at least 70%. Statistical 
analysis was performed by binning the time series at each temperature. Fig. IHl shows 
snapshots of two low-temperature configurations with A^ = 12 x 12 generated by the 
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Figure 6. Snapshots of configurations generated by the exchange method at T « 
iO~^ J on a 12 X 12 lattice for planar spins and J > 0. Periodic boundary conditions 
are imposed at the edges. 

exchange method with planar spins. In contrast to the case J < 0, it is difiicuh to 
decide on the basis of such snapshots if any order arises for J > 0: there are definitely 
large fluctuations including domain walls in the system. These fluctuations are in fact a 
necessary ingredient of the order-by-disorder mechanism. A careful quantitative analysis 
is therefore clearly needed. 

We start with the specific heat C, shown in Fig. There is a broad maximum at 
T ^ 0.25 J for planar spins and T ps; 0.3 J for spherical spins (see Fig. [7|^a)). However, 
this maximum does not correspond to a phase transition, as one can infer from the 
small finite-size effects. There is a second small peak in C at a lower temperature 
T ^ 0.02 J, see panels (b) and (c) of Fig. [3 Since this peak increases with A^, it is 
consistent with a phase transition. Note that at the lowest temperatures C/N is clearly 
smaller than 1/2 and 1 for planar and spherical spins, respectively. Indeed, in-plane 
fluctuations around the 120° state yield one branch of soft modes, which we expect 
to contribute only iV/12 to C rather than as the other two branches. Thus, we 

expect C/N = 5/12 = 0.41666 ... for planar spins and C/N = 11/12 = 0.91666 ... for 
spherical spins in the limit T ^ 0. Our MC results tend in this direction, but we have 
not reached sufficiently low temperatures to fully verify this prediction. Lastly, we note 
that the difference between C/N for planar and spherical spins is consistent with 1/2 
within error bars only for T < 0.02 J. 

Next, in Fig. IHlwe show the square of the sublattice magnetization (jHl) for J > 0. 
First, we observe that it remains small for almost all temperatures and shoots up just 
at the left side of Fig. |H1 where we measure a maximal value of ~ 0.4 for T = J/100 
on the = 6 X 6 lattice. This is consistent with a phase transition into a three- 
sublattice ordered state in the vicinity of the low-temperature peak in the specific heat 
C. In marked difference with the case J < (see Fig. |3]), here we observe only small 
differences between planar and spherical spins (at least for the sizes where we have data 
for planar spins, i.e., A^ = 6 x 6, 9 x 9 and 12 x 12). 



Finite-temperature ordering in a two-dimensional highly frustrated spin model 8 




T/J T/J 



Figure 7. Exchange MC results for the specific heat C for J > 0. (a) Full lines are for 
planar spins; dashed lines for spherical spins. Increasing line widths denote increasing 
system sizes N. Error bars are on the order of the width of the lines in this panel. The 
other panels show the low-temperature parts of the results for spherical spins (b), and 
planar spins (c). The system sizes in panels (b) and (c) range from = 6 x 6 (bottom) 
to = 18 X 18 (panel (b), top), and = 12 x 12 (panel (c), top), respectively. 

Finally, Fig. El shows our results for the Binder cumulant, as defined in (0}. In 
contrast to the sublattice magnetization m^, the Binder cumulants for planar and 
spherical spins are close to each other only for T < 0.1 J, provided we scale to the 
same prefactor A in the definition Again, our best estimate for the transition 

temperature Tc is obtained from the crossings of the Binder cumulants at different sizes 
N shown for planar spins in the inset of Fig. IHl The only meaningful crossings are those 
of the iV = 6 X 6 with the = 9 x 9 and 12 x 12 curves, respectively. This leads to a 
rough estimate 

yplanar ^ Q Q^25 J. (6) 

This estimate is indistinguishable from the corresponding value for spherical spins [14], 
as is expected since the spins lie essentially all in the x-y-p\ane for such low temperatures. 




Figure 8. Exchange MC results for the square of the sublattice magnetization for 
J > 0. Full lines are for planar spins; dashed lines for spherical spins. Increasing line 
widths denote increasing N. Error bars are negligible on the scale of the figure. 
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Figure 9. Exchange MC results for the Binder cumulant for J > 0. (a) Full lines 
are for planar spins; dashed lines for spherical spins. Increasing line widths denote 
increasing N. Error bars are at most on the order of the width of the lines, (b) Binder 
cumulant for planar spins close to the critical temperature (jSJ. Lines of increasing 
slope are for increasing system sizes A^ = 6x6,9x9, and 12 x 12. 



5. Conclusions and outlook 

In this paper, we have investigated finite-temperature properties of the classical version 
of the Hamiltonian on a triangular lattice. 

For J < 0, we have clear evidence for a low-temperature phase with 120° Neel 
order. The transition temperature is of order \J\ and can be determined with reasonable 
accuracy. Planar and spherical spins yield qualitatively similar results, but there 
are quantitative differences, in particular the transition temperature Q is higher for 
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planar spins. We have performed further simulations for spherical spins [14] in order 
to determine critical properties. On the one hand, so far we have no evidence for a 
discontinuity at Tc, on the other hand the assumption of a continuous phase transition 
yields very unusual critical exponents [14] which violate the hyperscaling relation. Thus, 
the most plausible scenario may be a weakly first-order transition. 

For J > 0, the groundstates are macroscopically degenerate [14]. A thermal 
order-by- disorder mechanism [18] predicts the selection of another 120° ordered state. 
The corresponding phase transition is just at the limits of detectability even with the 
exchange MC method [19,20]: we find an extremely low transition temperature © 
which is two orders of magnitude smaller than the overall energy scale J. In this case, 
the spins lie essentially in the x-y-plane in the relevant temperature region such that 
we obtain quantitatively extremely close results for planar and spherical spins in the 
vicinity of Tc, apart from a constant offset in the specific heat. 

The features observed e.g. in the specific heat for the classical variant resemble 
those found in the original quantum model [12-14]. Since only much smaller systems 
are numerically accessible for the quantum model, the results for the classical variant 
are an important tool for understanding the quantum case. In particular, the finite- 
temperature phase transitions should be universal and may therefore be characterized 
in the classical model. However, highly accurate data is needed for that purpose. Our 
result that spherical and planar spins exhibit the same qualitative features for J < 0, 
and that for J > the z-component is even quantitatively very small in the relevant 
temperature range, may be useful in this context. Namely, one may choose planar spins 
for further simulations, thus reducing the number of degrees of freedom to be updated. 
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